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Abstract 

We study bottomonium spectral functions in the quark-gluon plasma 
in the T and % channels, using lattice QCD simulations with two flavours 
of light quark on highly anisotropic lattices. The bottom quark is treated 
with nonrelativistic QCD (NRQCD). In the temperature range we con- 
sider, 0.42 < T/Tc < 2.09, we find that the ground states survive, whereas 
the excited states are suppressed as the temperature is increased. The po- 
sition and width of the ground states are compared to analytical effective 
field theory (EFT) predictions. Systematic uncertainties of the maximum 
entropy method (MEM), used to construct the spectral functions, are dis- 
cussed in some detail. 
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1 Introduction 



Quarkonia (heavy quark-antiquark bound states) are among the most important 
probes of the hot medium created in relativistic heavy ion colhsions. Unhke 
hght quarks, heavy quarks are predominantly created in the primordial hard 
collisions, and do not reach chemical equilibrium with the medium. Since J/ip 
suppression was proposed in 1986 cLS cl SI gnature of the formation of the quark- 
gluon plasma pQ, the charmonium system has been investigated intensively, both 
experimentally [2113] and theoretically . With the advent of the Large Hadron 
Collider, there has been increasing interest in bottomonium states as well, since h 
quarks are now for the first time being produced copiously in heavy ion collisions. 
In particular, the recent results from CMS indicate the survival of the T(15') 
state, but suppression of the T(25' + ?>S) states [6] (see Ref. |7] for results from 
STAR). A number of phenomenological studies have since attempted to explain 
this suppression pattern [8l[9] . It is generally expected that bottomonium provides 
a cleaner probe than charmonium, since statistical recombination of independent 
quarks and antiquarks plays a much less important role, and also since h quarks 
retain their nature as 'hard probes' to a larger extent. Cold nuclear matter effects 
are also expected to be simpler for h than for c quarks. 

Theoretically, quarkonium suppression has traditionally been investigated 
with potential models (see e.g. Refs. [101111] and references therein) and with lat- 
tice QCD computations of quarkonium spectral functions [T2l - [T9] . Exploiting the 
strongly-coupled nature of the quark-gluon plasma, studies using gauge-gravity 
duality have also been used, see e.g. Refs. [201 - I22] for recent results. Other studies 
can be found in Refs. [23|[2^ 

In the past few years, the theoretical understanding of quarkonium melting 
and in-medium modification has been improved substantially by casting the prob- 
lem in the language of effective field theories (EFTs) [251435] . By relying on scale 
separation between the heavy quark mass M and the temperature T of the quark- 
gluon plasma and on weak coupling to distinguish the inverse system size Ma^, 
the binding energy Mai ^"^^ ^^e inverse Debye screening length rriD y/a^T, a 
series of EFTs can be written down. One of the main outcomes of this formula- 
tion is the appearance of a complex heavy quark potential, where the imaginary 
part is generated by integrating out thermal fluctuations. For complex potential 
model studies, see e.g. Refs. [3S1EZ] as well as those listed above. Attempts at 
extracting the complex potential from lattice QCD can be found in Refs. [551155]. 

Various sequences of EFTs can be constructed, depending on the ordering of 
the scales. For instance, Refs. [251127] have focused on high temperature, where 
the bound state is about to fall apart, employing 



The corresponding bound state spectral functions are then considerably affected 
by the presence of the quark-gluon plasma. On the other hand, in Ref. [31] lower 



M > T > Mas >mD'> Ma]. 
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temperatures are considered, using 

M > Mas > T > Ma] > niD- (1.2) 

In this case the ground states are less affected and thermal effects can be cast in 
terms of thermal mass shifts and widths. 

In all cases, integrating out the heavy quark mass scale M yields nonrelativis- 
tic QCD (NRQCD) as an effective field theory. The appearance of further EFTs 
depends on the ordering of the scales and weak coupling. Since the applicability 
of weak coupling arguments is not guaranteed for temperatures up to 2 — ST^, 
where is the transition temperature between the hadronic phase and the quark- 
gluon plasma, it would be desirable to solve NRQCD in the quark-gluon plasma 
nonperturbatively, using lattice QCD simulations at finite temperature. This 
programme was recently initiated by usQ 

In Ref. [H] we studied S and P wave bottomonium correlators at four differ- 
ent temperatures, T/T^ = 0.42, 1.05, 1.40 and 2.09, using NRQCD for the heavy 
quark dynamics and relativistic two-flavour lattice simulations for the quark- 
gluon system. The main result of that analysis was the presence of strong tem- 
perature dependence in the P wave correlators (in the Xbi,b2,b3 channels), indicat- 
ing a melting of P wave bound states in the quark-gluon plasma. At the highest 
temperature, we found the behaviour of the P wave correlators to be consistent 
with nearly-free dynamics of the heavy quarks. 

The temperature dependence in the S wave correlators was much less vis- 
ible. The goal of this paper is to analyse in detail the 5" wave correlators, in 
the vector (T) and pseudoscalar (?7f,) channels, and construct the corresponding 
spectral functions at several temperatures in the hadronic phase and the quark- 
gluon plasma. Our main results can be seen in Fig. HJ which shows that as the 
temperature is increased the ground state peaks of the T and rji, remain visible, 
even though they broaden and reduce in height, while their excited states become 
suppressed at higher temperature and are no longer discernible at T/Tc ~ 1.68. 
The temperature dependence of the position and width of the ground state peaks 
is compared to analytical predictions obtained within the EFT formalism |31j . 
We note that the survival of the T(IS') state and suppression of excited states is 
consistent with the recent CMS and STAR results [HlCl. 

This paper is organised as follows. In the following section we decribe NRQCD 
as an effective field theory for QCD, focusing on finite temperature aspects. We 
discuss in particular how the nonrelativistic formulation has several advantages 
compared to standard relativistic dynamics at nonzero temperature. Lattice de- 
tails are collected in Sec. El The main part of the paper starts in Sec. HI where 
the high-precision euclidean correlators in the T and rjh channels are presented. 
The corresponding spectral functions are shown in Sec. [5l Here we also compare 

^See Ref. [3D] for an early pioneering study. 
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our results with analytical EFT predictions. A discussion of the maximum en- 
tropy method |12], used to construct the spectral functions, and of systematic 
uncertainties is given in Sec. [6l We summarise in Sec. [71 The Appendix contains 
an analysis of lattice artefacts in NRQCD spectral functions in the absence of in- 
teractions. Some preliminary results have previously been presented in Ref . |13] . 

2 NRQCD at nonzero temperature 

NRQCD is an effective theory of QCD where physics above the scale of the 
heavy quark mass is integrated out [^H - ffFj . It differs from heavy quark effective 
theory (HQET) in that terms in the NRQCD lagrangian are ordered in powers 
of f = |p|/M, the typical velocity of a heavy quark in the heavy quarkonium 
rest frame. In principle, there are infinitely many terms in such an expansion 
and taking into consideration more terms would mean more accurate relativistic 
corrections. However, in practice, only a small number of terms is necessary 
for a given accuracy, since v"^ is small (~ 0.1 for the bottom quark) and the 
series converges quickly. Also, including more terms in an effective theory would 
normally mean tuning more coefficients, resulting in a loss of predictive power. 
In practice, the coefficients of the NRQCD lagrangian can be calculated using 
perturbation theory since M is large (~ 5 GeV for the bottom quark) and often 
tree level values are sufficient when simulation parameters are chosen judiciously. 

In this work, we use the following 0{v^) euclidean NRQCD lagrangian density 
for the bottom quark [46] . 

£ = £o + 5A (2.1) 

with 
and 

+C2^ [V^Ud-e-e-d)v^ + x^(d-e-e-d)x] 

-C3^ [^^CT ■ (D X E - E X D) + ■ (D X E - E X D) x] 
-C4^[V^V-B7A-xV-Bx]. (2.3) 

Here Dr and D are gauge covariant temporal and spatial derivatives, ip is the 
heavy quark and x is the heavy anti-quark. The coefficients q = 1 at tree level. 

In contrast to the relativistic theory, the time evolution for nonrelativistic 
heavy quarks is an initial value problem. In particular, the presence of a nonzero 
temperature is not imposed as a boundary condition in the temporal direction 
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of the heavy quark field. Instead, the effects of temperature enter for the heavy 
quarks as they propagate through the thermal medium of light quarks and gluons. 
Since we are considering temperatures T <^ M, we expect that finite temperature 
can be taken into account without affecting the effective field theory nature of 
NRQCD. 

The absence of thermal boundary conditions simplifies spectral relations con- 
siderably. In the relativistic formulation, the euclidean correlator and its spectral 
function are related via 

G{t) = ^K{t,co)p{u;), (2.4) 

with the kernel 

cosh[a;(r-l/2T)] 

sinh(a;/2T) ' ^^.S) 

Temperature dependence enters in two ways: kinematically due to the periodic 
boundary conditions, refiected in the periodicity of the kernel, and dynamically 
due to the propagation through a temperature-dependent medium. It is impor- 
tant to disentangle these two, since the first one is present even in the absence of 
interactions and does not refiect the effects of the thermal medium. 

In NRQCD the kinematical temperature dependence is absent. This can be 
seen in a number of ways. Following Ref . [27] , we write u = 2M + u' and drop 
terms that are exponentially suppressed when M ^ T. The spectral relation 
(12. 4p then reduces to 

— exp{-u'r)p{u') (NRQCD), (2.6) 

even at nonzero temperature. This reflects the fact that the NRQCD propagator 
is constructed from an initial-value problem. Physically it implies that the heavy 
quarks are not in thermal equilibrium with the light-quark-gluon system, but 
instead appear as probes. 

This simpliflcation also removes the problems associated with the so-called 
constant contribution |l8]. In the small energy limit, the product of the relativis- 
tic kernel and the spectral density is independent of euclidean time 



lim i^(r, u;)p(a;) =2r^ , (2.7) 



uj=0 



where we used that the spectral function p{uj) is an odd function in u and in- 
creases linearly at small u. This is relevant for transport coefficients [H] and 
for conserved charges, in the presence of which spectral functions will have a 
contribution of the form 

p{uj) = x27r(X'5(u;) + contribution at larger u, (2.8) 
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where x is a susceptibility. Spectral weight at vanishing energy will therefore 
yield a constant, additive contribution to the euclidean correlator. It has been 
argued [48l[50] that this constant contribution interferes with the interpretation 
of charmonium survival or melting, as seen by lattice QCD simulations [T214T9] . 
It also requires a modification [51] of Bryan's algorithm in the implementation 
of the maximum entropy method |l2l[52]. In NRQCD, the contribution at small 
energies is absent, since only energies above 2M are presently In summary, in 
the heavy quark limit the spectral relation simplifies considerably, removing the 
problems associated with thermal boundary conditions. All temperature effects 
seen in the correlators are thus due to changes in the light-quark-gluon system. 

3 Lattice formulation 

We solve NRQCD nonpertubatively using lattice QCD, and let the bottom quarks 
propagate through a medium of gluons and two flavours of light quark. Gauge 
configurations with Nf = 2 dynamical light Wilson-type quark flavours are pro- 
duced on highly anisotropic lattices [C, = ag/ar = 6) of size A'^^ x N^. A summary 
of the lattice datasets is given in Table [H while more details of the lattice ac- 
tion and parameters can be found in Refs. [T6l[55]. In the light quark sector, 
m.j^/mp ~ 0.54, which implies that the light quark masses are comparable to the 
strange quark mass. 





Nr 


T(MeV) 






12 


80 


90 


0.42 


250 


12 


32 


230 


1.05 


1000 


12 


28 


263 


1.20 


1000 


12 


24 


306 


1.40 


500 


12 


20 


368 


1.68 


1000 


12 


18 


408 


1.86 


1000 


12 


16 


458 


2.09 


1000 



Table 1: Summary of the lattice data set. The lattice spacing is set using the IP— 
IS spin-averaged splitting in charmonium [17], corresponding to ~ 0.162fm, 
~ 7.35 GeV. The anisotropy is as/a^ = 6. 

In Ref. [H] we studied the effect of temperature on S and P wave bottomo- 
nium correlators. Compared to that study, we have greatly increased the number 
of configurations and expanded the number of temperature values (in Ref. [H] 
we only considered T/Tc = 0.42, 1.05, 1.40 and 2.09). We have also improved the 
nonperturbative tuning of the bare anisotropy in the actionlf] 

^ Heavy quark diffusion has been studied using heavy quark effective theory [551IM] . 
All data analysed here correspond to Run 7 in the terminology of Ref. [16] . 
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In order for NRQCD to be a consistent effective field theory in a lattice sim- 
ulation, the lattice spacing a^, acting as a short-distance cutoff, has to be kept 
finite and satisfy Mas ~ 1- Finite lattice spacing errors can then be systemat- 
ically improved as they appear at the same order as relativistic effects due to 
Qs ~ 1/M. There are many equivalent discretizations of the continuum NRQCD 
lagrangian density discussed above. Following earlier studies of heavy quarko- 
nium spectroscopy at zero temperature [56] - [58] . we calculate the heavy quark 
Green function on an anisotropic lattice using 

G(x,r = 0) =5(x), 
G(x, r = a.) = (^1 - f/] (x, 0) (l " ^) " ^'(x, 0), 

G(x, r + a.) = (^1 - f/i(x, r) (^l - - 5H) G(x, r), (3.1) 

where 5'(x) is the source, the lowest-order hamiltonian reads 

A(2) 

Ho = , (3.2) 

° 2M' ^ ^ 



and 



6H = ^ + ^(A . E - E ■ A) - -^a- ■ (A x E - E x A) 

g ^ a2A(4) a,(A(2))2 
"'^•B + ^TTTT T^-JT^- (3-3) 



2M 24M 16nM2 

The integer n controls the high-momentum behaviour of the evolution equation. 
Since the bottom quark is heavy enough, we take n = 1. The last two terms in 
6H are corrections to the kinetic energy term at finite lattice spacing |15]. The 
lattice covariant derivatives are defined as 



2as 

A(2)^ = Af V = J2~^ ^Ui{x)'^{x + i) - 2'^{x) + f//(x - t)ilj{x - i) 
A(^)^ = $^(Af ))^^, (3.4) 

i 

and E and B in Eq. (13. 3p are lattice definitions of the chromoelectric and chro- 
momagnetic fields. We use tadpole improvement [59] : 

U,{x) ^ Uo{x) ^ (3.5) 

Us Ur 
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where 

lis T cirG the average space-hke (s) and time-hke (r) hnks, determined from 
the plaquette expectation values; although in practice, the time-like mean link 
Ur is set to 1. The coefficients q are then set to 1. Note that Ug ^ Ur, since the 
lattice is anisotropic. 

An accurate determination of bottomonium spectroscopy requires careful tun- 
ing of the bare heavy quark mass M to satisfy NRQCD dispersion relations [56] . 
To study the finite-temperature modification of NRQCD propagators, an approx- 
imate choice of a^M such that M ~ 5 GeV is sufficient. 

There are many sources of systematic error in a lattice NRQCD calculation 
|58] . The three usual ones are effects of finite lattice spacing, finite volume 
and light quark vacuum polarization. Two more arise from the effective field 
theory nature: relativistic effects and radiative corrections to the couplings in the 
NRQCD expansion. Among these, finite volume effects are expected to be small 
for bottomonium since the physical size of bottomonium is small. Relativistic 
effects from the neglected higher order terms beyond 0{v^) are expected to be 
small as well since f ^ ~ 0.1. From experience in heavy quarkonium spectroscopy 
[56|[57]. tadpole improvement is expected to reduce radiative corrections. The 
light quark masses used in this work are somewhat large but since the effects of the 
light quark vacuum polarization in heavy quarkonium physics is small, systematic 
effects from this are expected to be minor. The presence of a finite lattice spacing 
with the condition Mag ~ 1 is a well-known issue in lattice NRQCD, ruling out 
a continuum limit of NRQCD results, but at finite temperature this problem is 
no worse than at zero temperature. In all, the qualitative features of the finite 
temperature behaviour of bottomonium reported in this work are expected to 
remain valid even after the careful consideration of systematic errors. 

4 Correlators 

The starting point for the remainder of this paper is the high-precision euclidean 
correlators in the vector (T) and the pseudoscalar (r/^) channel, obtained by solv- 
ing the NRQCD evolution equations. In Fig. [1] (top) we show these correlators, 
normalised with the value at r = and on a logarithmic scale, for all tempera- 
tures available. Due to the use of NRQCD, there is no periodicity in the temporal 
direction. Combined with the large anisotropy, this implies that many temporal 
lattice points are available for the analysis, even at the highest temperature. 

It is clear from the plots on the top that the temperature dependence is 
very mild. In order to make the temperature dependence visible, we show on 
the bottom the ratio of the high-temperature correlators with the one in the 
hadronic phase (T/T^ = 0.42, A^^^ = 80). We observe that the effect of increasing 
the temperature is monotonic and always below 3%. We remark here that the 
ratios depend on the sources used in the NRQCD evolution; the results shown 
here are obtained with point sources. Since potential model calculations typically 
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5| (vector) 
Upsilon 



16 



32 



A r/r=o.42 

c 

• r/r=i.05 

< TIT -\. 20 




48 



64 




16 



32 



48 



64 



80 




Figure 1: Euclidean correlation functions G{t) as a function of the euclidean 
time in the vector (T) channel (left) and the pseudoscalar (77^) channel (right), 
for all temperatures available, using point sources. At the top the correlators are 
normalised with the value at r = and shown on a logarithmic scale, while on 
the bottom the high-temperature correlators are normalised with the correlator 
at the lowest temperature, T/Tc = 0.42 (A',- = 80). The errors are smaller than 
the symbols. 

use point (delta function) sources as well, a comparison between potential model 
predictions and our NRQCD results should be possible. 

In contrast to the relativistic case, the temperature dependence seen here does 
not arise from the thermal boundary conditions, but solely from the presence of 
the medium of gluons and light quarks at different temperatures. In terms of the 
spectral relation, 

G{t) = j ^e-Xo;), (4.1) 

this is reflected in the temperature-independent kernel e"'^'^. In the relativistic 
case, the temperature dependence of the kernel means a direct comparison be- 
tween correlators at different temperatures is not straightforward. Often the kine- 
matical temperature dependence is eliminated by using so-called reconstructed 
correlators, which requires the calculation of a spectral function at the lowest 



9 




Figure 2: Relative error err [G(r)]/G'(r) as a function of euclidean time in the 
vector (left) and the pseudoscalar (right) channels. The data at A^,- = 80 and 24 
are rescaled to take into account the lower number of configurations. 

available temperature ^U\. This is not needed here and therefore the ratios 
in Fig. [T] (below) are a proper reflection of the presence of the temperature- 
dependent medium. 

The statistical errors are small and not visible in the plots discussed above. 
To illustrate this, we show the statistical relative error, i.e., err[G(r)]/G(r), as a 
function of euclidean time in Fig. |2l At all but two temperatures, there are 1000 
configurations available. To take this into account, we rescaled the = 80 data 
by a factor 2 (250 configurations) and the A^,- = 24 data by a factor \/2 (500 
configurations). We then observe that the relative error is of the order of 10~^ 
and that it increases as the temperature of the quark-gluon plasma is increased, 
indicating larger thermal fluctuations in the hot phase. 

From the data at the lowest temperature, we extract the masses of the ground 
states and the first excited states using standard exponential fits. The results are 
summarised in Table |2] and were already presented in Ref. HI] El In NRQCD, 
all energies are determined only up to an additive constant, E = Eq + AE. 
Taking the T(IS') mass from the Particle Data Book [ni] to set the scale, we find 
Eo = 8.57 GeV. 

"'in Ref. 01] a temporal lattice spacing of = 7.23 GeV rather than 7.35 GeV was used, 
correcting this results in a small change in the mass predictions in the third column. 
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state 




Mass (MeVj 


bxp. (MeVj IdJJ 


^ So{r]b) 


0.118(1) 


9438(8) 


9390.9(2.8) 


2^So{Vb{2S)) 


0.197(2) 


10019(15) 






0.121(1) 


9460* 


9460.30(26) 


23^1 (T') 


0.198(2) 


10026(15) 


10023.26(31) 




0.178(2) 


9879(15) 


9898.3±1.1+};? [62] 




0.175(4) 


9857(29) 


9859.44(42)(31) 




0.176(3) 


9864(22) 


9892.78(26)(31) 




0.182(3) 


9908(22) 


9912.21(26)(31) 



Table 2: Zero temperature bottomonium spectroscopy. The l'^S'i(T) state is used 
to set the scale. Here we concentrate on the T and r/^ states. 

5 Spectral functions 

We extract spectral functions from the euclidean correlators presented above 
using the Maximum Entropy Method (MEM) A straightforward inversion 
of Eq. (14. ip is not possible, since euclidean correlators are determined numerically 
at a finite number of points, whereas spectral functions are in principle continuous 
functions of the energy u. Using the ideas of Bayesian probability theory, one 
can construct the most probable spectral function by maximizing the conditional 
probability P[p\DH], where D indicates the data and H some additional prior 
knowledge, encoded in a default model. In this section we present the results, 
while a discussion of the systematic uncertainties is given in the next section. We 
only consider spectral functions at zero spatial momentum. 

The results for the spectral functions at the lowest temperature are given 
in Fig. 131 The vertical lines indicate the position of the ground state and first 
excited state from Table El We observe that the ground state appears as a very 
narrow peak in the spectral function, while the first excited state is broader 
and overlaps with structure at larger energy. The third feature, visible in the 
inset, is presumably a combination of higher excited states and lattice artefacts 
(see Appendix |A])|3 We note that quadruple precision is required for the MEM 
inversion, due to the large number of points (A'^ = 80) and the exponential 
fall-off. 

The main result of this paper is the spectral functions at the various temper- 
atures, shown in Fig. |4] for the vector (T) channel (upper panel) and the pseu- 
doscalar (77^) channel (lower panel). Every panel contains two adjacent tempera- 
tures, from the coldest (T/Tc = 0.42) in the top left to the hottest (T/Tc = 2.09) 
in the bottom right. In each panel, the lower temperature is depicted with a 

^Note in particular that it cannot be identified with the second excited state T{3S), with a 
mass of 10.3552 GeV [SI]. 
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Figure 3: Spectral functions p{u}), normalised with the heavy quark mass, as 
a function of energy at the lowest temperature in the vector (T) and the pseu- 
doscalar (t]^) channels. The vertical lines indicate the positions of the ground 
and first excited state obtained via standard exponential fits. The insets show a 
close-up. 



dashed line and the higher one with a solid line. As the temperature is increased, 
we observe that the ground state peak remains visible, even though it broad- 
ens and reduces in height. The excited states become less pronounced as the 
temperature increases and are no longer discernible as a separate peak between 
1.4 < T/Tc < l.Gsj^ This can be interpreted as the survival of the 15" states, 
but a melting or suppression of the excited states. From the analysis in Sec. |6l 
we note here that we consider the results for the spectral functions to be robust 
for all temperatures, with the exception of the two highest temperatures where 
uncertainties due to the limited statistics and euclidean time range remain. 

We note that the area under the curve is determined by the sourc dll at r = 
and the spectral relation, 

J ^p{tu,p = 0) = J d^xG{r = 0,^) = Jd^xS{x.), (5.1) 

and is independent of the temperature. 

The peak position E and width F of the ground states can be extracted by 
fitting the peaks to a Gaussian function. We fit to the left side of the peak, 
to avoid contamination from the features at larger u. In Fig. |S] (top panel) 
we show the temperature dependence of the mass shift AE, normalised with 
the heavy quark mass. Recall that in NRQCD only energy differences can be 
determined, and that the total energy is E = Eq + AE, where Eq = 8.57 GeV in 

^We also note that the second peak immediately above Tc is presumably a combination of 
the first excited state and other features. 

^The point source is defined to be unity for each of the upper two-component spinor indices 
as well as for each of the colour indices. 
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CO [GeV] 




9 10 11 12 13 14 9 10 11 12 13 14 9 10 11 12 13 14 15 



CO [GeV] 

Figure 4: Spectral functions p(w), normalised with the heavy quark mass, in 
the vector (T) channel (upper panel) and in the pseudoscalar (77^) channel (lower 
panel) for all temperature available. The subpanels are ordered from cold (top 
left) to hot (bottom right). Every subpanel contains two adjacent temperatures 
to facilitate the comparison. 

our case. The temperature dependence of the width is shown in Fig. [5] (bottom 
panel). Note that the width is normalised with the temperature. The error bars 
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0.19 



0.19 r 




Figure 5: Position of the ground state peak AE, normalised with the heavy 
quark mass (upper panels), and the upper limit on the width of the ground state 
peak, normalised with the temperature (lower panels), as a function of T/T^ 
in the vector (T) and the pseudoscalar (77^) channels. The error bars denote 
the systematic uncertainty with the left error bars representing the error from 
the finiteness of the last time in the fitting window, r2, and the right error bar 
representing the error from the finite statistics (see Sec. [6]). The lines in the upper 
plots indicate expected analytical behaviour assuming weak coupling above T^. 



indicate systematic uncertainties in extracting the peak position and width from 
the peaked structure. In Sec. [6] these uncertainties are discussed in detail. Based 
on this discussion we conclude conservatively that the width shown in Fig. [5] is 
better interpreted as an upper bound, rather than the width itself. 

To see whether these results are reasonable, we now take them at face value 
and contrast them with analytic predictions derived assuming a weakly coupled 
plasma. According to Ref. [31], the thermal contribution to the width is given, 
at leading order in the weak coupling and large mass expansion, by 

^ = ^-^al ^ U.27al (5.2) 
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i.e., the width increases hnearly with the temperature |f| If we take as an estimate 
from our results that F/T ~ 1, we find that this corresponds to ~ 0.4, which is 
a reasonable result. In the same spirit the thermal mass shift is given in Ref. [31] 
by 

17 IT 

^^thermal = ~^^^J^ " 5.93a^ — . (5.3) 

In these simulations we have ~ 220 MeV, M ~ 5 GeV. Taking these values 
together with as ~ 0.4 as determined above, Eq. fl5.3p becomes 

In order to contrast our results with this analytical prediction, we have compared 
the temperature dependence of the peak positions to the simple expression 

_ = , + 0.0046 (-j . (5.5) 

where c is a free paramater. This is shown by the dashed lines in Fig. [5] (top 
panel). While we are not in a position to confirm or rule out the quadratic tem- 
perature dependence due to the systematic uncertainties in the MEM analysis 
and the fitting of the ground state peaks, we note that our results are not incon- 
sistent with this. In particular, the absolute scale of temperature variation seems 
to be of the correct order. We also note that the mass just above is reduced 
with respect to the low temperature one. 

6 Systematics and uncertainties 

In order to assess the robustness of the results presented above, in this section we 
discuss the main systematic uncertainties, namely the dependence on the default 
model, the choice of cjmin, the number of configurations and the euclidean time 
window, and explain the error estimates of the previous section. We only show 
results for the vector channel; the ones in the pseudoscalar channel are similar. 

Default model 

The MEM procedure includes a "default model" m(a;) which is used to define 
the entropy term. 



S 



2^ 



p{uj) — ■m{uj) — p{uj) In -^^^^ 



m(ui 



(6.1) 



*A linearly rising width with temperature was also predieted in Ref. 
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Figure 6: Dependence of the spectral function on the default model chosen in 
the MEM analysis for A^^ = 32 (left) and 16 (right), in the vector channel. The 
default models favoured by the minimization are denoted with a *. 

in the MEM approach. It is expected that with poor data the spectral function 
will resemble the default model, since this minimises the entropy, but that with 
precise data the choice of default model becomes irrelevant. The usual proce- 
dure is to choose a default model which has the same functional form as the 
(continuum) free spectral function in the channel under consideration. For the 
nonrelativistic S wave, this is m(a;) = moi/cJ; see Eq. flA.4p . An alternative is 



to use a constant default model, m{u) = mo. The constants mo can be set by 
minimising the between the data and the correlation function defined from 
the default model, 

Gdcf(r)= / —e~-^m{uj). 
I, , ■ ^vr 

To illustrate the absence of default model dependence in our analysis, we show 
here results for the following six default models: 

• m{uj) = mo-y/cJ, with mo/mQ = 0.1, 1, 10, 

• m(cj) = mo, with mo/mg = 0.1, 1, 10, 

where in both cases mg is determined by minimizing described above. 

In Fig.[6]we show the resulting six spectral functions, as well as the six default 
models used, at T/Tc = 1.05 and 2.09. We observe that the spectral functions 
do not resemble the default models and that there is almost no variation in the 
spectral functions, even with default models varying over two orders of magnitude. 
We conclude that there is no significant default model dependence in this analysis, 
even at the highest temperature. 
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TV 


TV 


T (vector channel) 


rjb (pseudoscalar channel) 


80 


4000 


0.00 


2.00 


0.00 


2.00 


32 


1000 


0.00 


1.50 


0.00 


2.00 


28 


1000 


-0.04 


1.46 


-0.08 


1.92 


24 


1000 


-0.10 


1.40 


-0.10 


1.90 


20 


1000 


-0.10 


1.40 


-0.10 


1.90 


18 


1000 


-0.10 


1.40 


-0.10 


1.90 


16 


1000 


-0.12 


0.88 


-0.12 


1.88 



Table 3: Details of the parameters used in the MEM analysis. Note that Au = 
(i^max — i^min) /N^^, and that ttrOo = and 2 correspond to 8.57 GeV and 23.3 
GeV respectively. 



Energy window 

Since in NRQCD the heavy quark mass scale is removed, the additive normaliza- 
tion of the energy scale has to be determined by comparing to a physical state. 
For this reason, setting the lower boundary Wmin = in the spectral relation 
( 14. ip is not a priori justified and it might be necessary to allow for the possibil- 
ity of negative energies (wmin < 0) in the energy integral. Indeed, in the high 
temperature cases we noticed that when Wmin ~ 0, the MEM spectral functions 
have a spike in the lowest-energy bin, which we interpret as an unphysical effect. 
Reducing Wmin to negative values reduces this spike until it is essentially absent 
when arOJmin ~ —0.10. At the lower temperatures (A'',- = 80,32), no spikes were 
observed with Wmin = 0. In Table [3] we list the parameters used in the MEM anal- 
ysis, denotes the number of points in which the energy interval 
is divided. 

Number of configurations 

In Fig. [7]we illustrate how the position and width of the ground state peak depend 
on the number of configurations used in the MEM analysis. For clarity, we only 
show results for a selection of Nr values. At each temperature, we divide the total 
number of configurations (1000) into 20 groups of 50, 10 groups of 100, 5 groups 
of 200, 3 groups of 333, and 2 groups of 500, and compute the spectral function 
for each of them. The resulting peak position and width are shown as a function 
of 1 / ^yNcigs- At the lower temperatures (A^- = 32,28,24) we observe that both 
are essentially independent of the number of configurations used and the results 
are stable. At the higher temperature (A^,- = 20) there is a larger spread for the 
low-statistics results, but the average results are again quite stable, with at most 
only a slight decrease in the width and mass as the statistics are increased. It is 
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Figure 7: Dependence of the position (AE) and width (F) of the ground state 
peak extracted from the spectral function on the number A^cfgs of configurations 
used, in the vector channel. 

only at the highest temperatures (A',- = 18, 16) that there is a clear decrease in 
the position and width with increased statistics. We conclude that the statistical 
error due to the finiteness of the ensemble does not prevent us from determining 
the ground state mass and width, at all but the highest two temperatures. 

To translate the dependence on the number of configurations into an uncer- 
tainty in the position and width of the ground state, we proceed as follows. The 
central value is taken from the analysis with 1000 configurations. The statistical 
error due to the finite number of configurations is defined by taking the difference 
between the central value and the value obtained with 500 configurations which is 
furthest from the central value. This statistical error is shown as the right-hand 
error bar in Fig. El 

Euclidean time window 

The euclidean time window included in the MEM analysis is ri < r < r2, with 
the initial time equal to the first time slice, i.e. ri = In Fig. [8] we show 

the dependence of the constructed spectral functions on T2 for A^,- = 20 (left) 
and Nt = 16 (right), varying T2 from a small value up to iV,- — 1. Both cases 
are at high temperature, T/T^ = 1.68 and 2.09 respectively. We observe that at 
A'^^ = 20 the result is stable and the limitation of having only a finite number of 
points in the temporal direction does not appear to be a problem. At the highest 
temperature, on the other hand, we note that the results have not yet stabilised, 
so that the uncertainty in constructing the spectral function is larger. 

In order to quantify this, we have determined the T2 dependence of the po- 
sition and peak of the ground state at all temperatures. The result is shown in 

^Varying ti within reason did not have a significant effect. 
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Figure 8: Dependence of spectral functions on the euclidean time window = 
Ti < r < r2, for several values of T2 for A^^ = 20 (left) and 16 (right), in the 
vector channel. 
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Figure 9: Dependence of the position {AE, left) and width (F, right) of the 
ground state peak extracted from the spectral function on the euclidean time 
window ttr = Ti < T < T2 used in the MEM analysis, in the vector channel. 

Fig. |9] for a selection of A',- values. For the lowest temperature {Nr = 80), we 
observe that AE is stable and consistent with the result from the standard fits 
directly to the euclidean correlator. Moreover, there is evidence that the width 
F decreases towards zero as T2 increases. The position and width at N^- = 32, 28, 
or T/Tc = 1.05, 1.20, appear to behave as in the hadronic phase, except for a 
possible flattening out of the width for Nr = 28 at large T2. A similar behaviour 
is found for A',- = 24, with a clearer tendency for both the position and width to 
flatten out at large T2- On the other hand, at = 20 both position and width 
are consistently above their low-temperature values, and may reach a plateau for 
large T2, consistent with Fig. [HI For N^- = 18 and 16, no stability is seen. For 
Nj. = 16, it follows from Fig. [S]that when T2 is too small, the ground state and 
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Figure 10: The probability -P(a) used in Bryan's approach for A^.^ = 32 (left) 
and 16 (right). The horizontal lines indicate the interval used in the integral in 
Eq. (ED. 



features at larger energies overlap and therefore the ground state peak cannot 
be resolved. At r2/a^ > 13, two structures become visible, which explains the 
rapid drop of the position and peak. We conclude therefore that at most tem- 
peratures above the position and width of the groundstate peak appear to be 
stable, except at the highest two temperatures where substantial dependence on 
T2 remains. 

To translate the finiteness of r2 into a systematic uncertainty on the position 
and width of the ground state, we adopt the following procedure. The central 
value is taken from the analysis with T2 = — 1. The systematic error due to 
a finite T2 is defined by taking the difference between the central value and the 
value obtained with T2 = — 2. This systematic error is shown as the left-hand 
error bar in Fig. [5l 



Bryan's approach 

In the MEM approach, Bayes' theorem implies that the entropy term S is bal- 
anced against the usual maximum likelihood x^. The function Q to be maximized 
is 

Q = aS-x\ (6.2) 

where 5* is defined in Eq. (16. ip and a is a normalisation constant which is to be 
determined. In Bryan's approach [52] the spectral function Pa{oj) is calculated 
for each value of a and the final spectral function is obtained by performing the 
convolution integral, 

Pfiniiiiuj) = j daP{a)pa{uj). (6.3) 
Here P{<y) is the probability that a is chosen correctly. 
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In Fig. [ini the probability -P(«) is plotted for the T channel for A',- = 32 and 
16. In practice, the limits in the integral in Eq. (16. 3 P are taken such that P{a) is 
greater than some fraction of its maximum value Pmax- In this case, we use the 
interval defined as -P(a)/Pmax > 0.1 and the resulting interval is shown in Fig. 
by the horizontal lines. In all cases we find that the probability is well defined 
with a clear maximum. 

7 Conclusion 

In this paper we analysed lattice QCD results for S wave bottomonium corre- 
lators (in the T and rjb channels). The heavy quarks are treated with NRQCD 
and propagate through a two-flavour quark-gluon medium at seven temperatures 
between 0.4Tc and 2.1Tc. Highly anisotropic lattices, with a^/ar = 6, are used to 
maximise the number of time slices available for the analysis. Spectral functions 
are constructed with the help of the maximum entropy method. 

The use of NRQCD has a number of advantages compared to relativistic 
quark dynamics. Since the presence of a nonzero temperature is not imposed as 
a thermal boundary condition, twice as many euclidean time points are available 
for the analysis compared to the relativistic case (at the same temperature and 
lattice spacing). More importantly, the spectral relation simplifies considerably, 
removing the problem with the so-called constant contribution in the correlator. 
Physically, it implies that all temperature dependence is due to the presence of 
the light-quark-gluon system. 

Our main results are spectral functions for the S waves, in the vector (T) and 
pseudoscalar (77^) channels. Our results suggest that the ground state survives 
up to the highest temperature we consider, whereas the excited states are sup- 
pressed and no longer visible at temperatures above lATc. We have extracted 
the position and width of the ground state peaks and found them to be consistent 
with analytical results obtained within the EFT framework, at leading order in 
the large mass expansion. 

Systematic uncertainties have been studied in some detail. For all tempera- 
tures, we found no dependence on the default model used in the MEM analysis. 
The position and width of the ground state peaks were shown to be stable as the 
number of configurations is increased, for all except the highest temperatures. 
Perhaps the biggest uncertainty comes from the finite number of euclidean time 
points, but our analysis suggests that the peak position and width can still be 
reliably determined for temperatures up to 1.7Tc. 

We hope that our results will be useful for further EFT and potential model 
studies. After constructing a temperature-dependent potential, one usually com- 
putes spectral functions and the corresponding euclidean correlators. It would 
be interesting to compare the outcome of such an exercise with our nonpertur- 
batively determined correlators. This would in particular be applicable to ratios 
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of finite-temperature correlators with the zero-temperature one, as in Fig. [H to 
cancel normalization factors and focus on the temperature dependence. 

Finally, we hope that our results will contribute to a further understanding 
of the recent experimental results for bottomonium in heavy ion collisions at the 
LHC and RHIC. 
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A Noninteracting lattice spectral functions 

In order to understand the effect of lattice artefacts, it is useful to construct 
lattice spectral functions in the absence of interactions, adapting the approach of 
Refs. [631E1] to lattice NRQCD. 

Let us start with free quarks in continuum NRQCD with energy Ep = p^/2M. 
The correlators for the S and P waves at zero spatial momentum are then of the 
form [27] 

G,M.2A,J^pV-.'^i%filV'\ (A,2) 



(27r)3 ^ 87r3/2 j 

where in the explicit evaluation we did not include an ultraviolet cutoff, since the 
integrals are finite for nonzero r. 
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This is easily expressed in terms of spectral densities, using 



G{t)= pM, (A.3) 

yielding 

ps{u) =47riV,y" ^^S{uj-2E^) = ^M^/V/^GM, (A.4) 

pp{co) = Attn, I -^y5{uj-2E^) = ^M'/'u'/'e{u). (A.5) 

Note that the minimal energy Wmin = corresponds to twice the heavy quark 
mass, due to the nonrelativistic approximation, \/p'^ + ^ M + Ep. In the 
presence of a momentum cutoff |p| < A, the maximum energy is finite and 
given by Umax = A^/M. Note also that there is no temperature dependence in 
the absence of interactions. All temperature effects enter via the propagation 
through the quark-gluon system. 

These results can easily be adapted to the lattice [63l|6l] , taking into account 
the lattice dispersion relation and the finite momentum integration over the first 
Brillouin zone. At lowest order in (unimproved) NRQCD, the dispersion relation 
is ^ ^ 

a.Ep = -logfl--^y (A.6) 



where M = QsM, ^ = ag/ar, and 

i=l 

with Ng the number of sites in a spatial direction. The lattice spectral functions 
then take the form 

Ps{u;)=^Y.^{uj-2Ep), (A.8) 
p 

p 

where the sums go over all momenta in the first Brillouin zone. Evaluating 
these numerically, as in Refs. |63l[6l], yields the spectral functions shown in 
Fig. [TTJ Here is taken large enough to be in the spatial thermodynamic 
limit, while there is no dependence on N^-. Improving the dispersion relation 
will give better agreement between lattice and continuum results at small u. 
The cusps result from reaching the edge of the Brillouin zone in the (1,0,0) or 
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Figure 11: Lattice spectral functions, in units of the temporal lattice spacing 
d-j- , as a function of Ct-o;, for S waves (left) and P waves (right) in anisotropic 
lattice NRQCD at lowest order, ignoring interactions, using ^ = ag/ar = 6 and 
aT-M = 0.75. The dashed lines indicate the continuum spectral functions in the 
absence of a cutoff. 

(1,1,0) direction (+ permutations). The maximal energy is determined by the 
maximal lattice momentum in the (1,1,1) direction, namely = 12, and equals 
a^ri^max — 0.503 for the parameters used here. Comparing these results with those 
for free relativistic quarks [63l|64], we conclude that the main difference is the 
temperature (or A^.^) independence. 
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